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Abstract. We present the results of a 2-D, two fluid (ions 
and neutrals) simulation of the ambipolar filamentation 
process, in which a magnetized, weakly ionized plasma is 
stirred by turbulence in the ambipolar frequency range. 
The higher turbulent velocity of the neutrals in the most 
ionized regions gives rise to a non-linear force driving them 
out of these regions, so that the initial ionization inhomo- 
geneities are strongly amplified. This effect, the ambipo- 
lar filamentation, causes the ions and the magnetic flux 
to condense and separate from the neutrals, resulting in a 
filamentary structure. 
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1. Introduction 

Magnetic fields contribute to the dynamical behavior of 
ionized astrophysical fluids such as those in the upper so- 
lar and stellar atmospheres, the interstellar medium and 
star-forming regions. Their influence is carried out by hy- 
dromagnetic waves which efficiently propagate perturba- 
tions, ensure a turbulent pressure or may even cause the 



development of instabilities (Arons & Max 1975) 



However, Kulsrud & Pearce (1969) showed that in the 
magnetized and weakly ionized interstellar medium hydro- 
magnetic waves are heavily damped in a frequency range 
(and thus scale) associated with ambipolar diffusion. At 
low frequency the neutrals are well coupled to the ions 
(which are tied to the magnetic field lines) and hydromag- 
netic waves propagate at the Alfven speed defined by the 
total inertia (given by ions+neutrals). At high frequency 
neutrals and ions are totally decoupled, and Alfven waves 
involve only the ions, which define a larger Alfven velocity. 
In the intermediate range (the 'ambipolar range', between 
the ion-neutral and neutral-ion collision frequencies z^ n 
and v n i) the neutrals are imperfectly coupled to the ions; 
this results in a drag which strongly damps the waves (see 



also Mclvor 1977, for a description of MHD turbulence in 
the partly ionized interstellar medium). 

The non-linear evolution of this process can cause 
an ambipolar filamentation of the magnetic field when 
a magnetized and weakly ionized plasma is stirred 
by hydromagnetic turbulence in the ambipolar range 
( Tagger et al. 1995 ) . If such a plasma presents small varia- 
tions in the ionization fraction (Z — pi/p n ), the turbulent 
velocity of the neutrals is higher in the most ionized re- 
gions, since they are better coupled to the ions. This gives 
rise to a force (given by the average of the v.Vv term) 
driving the neutrals out of the most ionized regions. By 
reaction the ions and the magnetic flux are compressed in 
these regions, so that the initial ionization inhomogeneities 
are strongly amplified. As a consequence a concentration 
of the flux tubes is expected to occur, producing a fil- 
amentary structure, so that turbulent energy would be 
converted into magnetic energy associated with the con- 
centration of the magnetic field. Tagger et al. (1995) pro- 
vided only order of magnitude estimates of the expected 
amplification of the ionization fraction. In this work we 
present a fully consistent 2-D non-linear numerical simu- 
lation of the mechanism in order to test its efficiency. 

The non-linear analysis is a fundamental tool to study 
the physics in certain astrophysical environments, such 
as molecular clouds, where the observed amplitudes of 
the turbulent velocities are comparable with the mean 
field velocities. The ambipolar filamentation mechanism 
might help to explain some well known problems aris- 
ing in magnetized, partially ionized astrophysical plas- 
mas. One of them is related with the observations of 
turbulence in molecular clouds. Observations show a fila- 
mentary structure, and strong supersonic motions result- 
ing in turbulent and magnetic energies in approximate 
equipartition, i.e., much larger than the thermal energy 
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( Myers fc Goodman 1988 ). The ambipolar filamentation 
mechanism would concentrate the magnetic field in in- 
tense flux ropes surrounded by essentially neutral clouds. 

Another possible application relates to the fibrilled 
structure observed in the magnetic field emerging from the 
solar photosphere, organized in very narrow flux tubes. 
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The ambipolar filamentation mechanism might provide 
an explanation for the spicules emerging from the photo- 
sphere: let us consider magnetic held lines raising from the 
photosphere. Then an Alfven wave of a given frequency, 
produced in the photosphere and initially below the local 
ambipolar frequency range, will propagate upward along 
the field lines and reach at high altitudes a plasma of much 
lower density, i.e., lower collision frequencies. It will thus 
be damped by ambipolar effects and can expel the neu- 
trals from the most ionized flux tubes, concentrating the 
magnetic flux in narrow tubes where strong vertical mo- 
tions can be expected. This would occur together w ith th e 
mechanism discussed by De Pontieu & Haerendel ( 1998|) . 
These prospects will be discussed in more detail in the last 
section of this work. 

We have carried out numerical simulations in which a 
weakly ionized and magnetized gas inside a cartesian box 
is submitted to a high amplitude oscillation emitted from 
one of its sides. The perturbation propagates inside the 
box as an Alfven wave with a frequency chosen to be in 
the ambipolar range, so that it will be strongly damped. In 
Section 2 we describe the dynamical equations that govern 
the evolution of a two fluid gas, together with the numeri- 
cal code and the boundary conditions used to solve them. 
We also discuss the numerical constraints present in our 
simulations. The results from the numerical experiments 
are presented in Section 3 and discussed in the context of 
the problems cited above in Section 4. 

2. The Numerical Code 

The magnetohydrodynamics (MHD) equations describing 
a two fluid (ions and neutrals) system are (Langer 1978): 



dpi 
dt 

dt 



+ V {pm) = 
+ V(p n v n ) = 



(1) 
(2) 



—^+(vi.V)vi = --g+{VxB)xB+p Pn {v n -v l )(i) 

dt pi 

— — + (v n .V)v n = g + PPi(Vi-v n ) (4) 

dt Pn 



dB 

~dt 



V x (vi x B) 



(5) 



For simplicity we assume an isothermal equation of 
state: 



Pi = PiC- s , 



Pn PnC s , 



(6) 
(7) 



where p, v and p are, respectively, the density, velocity 
and partial pressure of the ions (with subscript i) and neu- 
trals (with subscript n), g is the gravity, p is a constant 
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Fig. 1. Top: Initial ions (dotted line) and neutrals (solid 
line) density profiles in the z direction at y — 0. Bottom: 
Initial ions density profile in the y direction at z = 0. We 
use the grid zone indices to scale the horizontal (y) and 
vertical (z) directions in all figures. See Section 3 for an 
explanation of the parametrization used. 

such that Vi n = p,p n and v n i = ppi are the ion-neutral 
and neutral-ion collision frequencies, and c s is the sound 
velocity (assumed the same for ions and neutrals). We as- 
sume that ionization and recombination occur on a longer 
time scale than the one we consider. This should of course 
be checked for applications to specific astrophysical situ- 
ations. 

We have also checked that in these conditions the char- 
acteristics of the problems in which we are interested, 
namely the high electron densities in the case of solar 
spicules and the large spatial dimensions of molecular 
clouds, allow us to use these simplified two fluid MHD 
equations instead of the full set of the three fluids (elec- 
trons, ions and neutrals) equations which describe the dy- 
namics of a weakly ionized gas. 

In order to allow for the long time scales considered 
(hundreds of Alfven times), associated with numerical 
constraints giving a short time step, we simplify the prob- 
lem by making it two-dimensional: all quantities are x- 
invariant, and only the perturbed current has a component 
along x. Therefore all quantities depend only on z (the 
vertical direction) and y (the horizontal one) on a carte- 
sian grid. In this geometry the distinction between shear- 
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Alfven and magnetosonic waves disappears, but waves 
propagating along z retain the properties of the usual 
AlfVen waves, that is, they twist the field lines and to 
lowest order they are not compressional. 

The numerical code used in our simulations is based 
on the same general methods as the ZEUS-2D code 
pone fc No7 man 1992a|). It solves the MHD equations 



on a staggered mesh using the method of finite differ- 
ences with a time explicit, operator split scheme. Densities 
and pressures are zone centered quantities while the veloc- 
ity components are centered at their corresponding zone 
faces. As in ZEUS-2D the solution procedure is arranged in 
two main steps, first taking into account the source terms 
in the right-hand side of the equations, and then solving 
for the advection terms in the left-hand side. The main 
difference with ZEUS-2D is the treatment of the gravita- 
tional and magnetic terms. An initial state of hydrostatic 
equilibrium is assumed, by fixing the vertical density and 
pressure profiles, from which we derive a value for the 
gravity; this results in an ad-hoc gravity profile leaving us 
some freedom to optimize the density profiles, limiting the 
numerical difficulties associated with the emission of the 
wave (see below) . The 2-D geometry allows us to describe 
the magnetic field in cqs. (||) and (||), in terms of the flux 
function ip as follows: 



B = B e z + e x x Vijj 

that the induction equation, eq. (|^) becomes: 



so 

d_ 

dt 



-VyB 



(8) 



(9) 



The momentum and induction equations are thus first 
solved without the advection terms, using time explicit op- 
erator split schemes. The MOC-CT algorithms included in 
ZEUS-2D flStone fc Norman 1992b] ) are not used to evolve 
the magnetic terms, since they are not required by this 
simple problem. Tests were carried out to verify the cor- 
rect propagation of the waves. 

The second step in the numerical code, still following 
ZEUS-2D, solves finite difference versions of the integral 
equations coming from the advection terms in eqs. (Q), 
§, § and (H): 



d 
dt 



— pi dV = - p t Vi. dS 



dV 



d 
dt 



37 / Pn dV = 



p n v n . dS 



dV 



d 
dt 

d 

dt 



Pi VidV 



p n v n dV 



PiViVi. dS 



d\' 



dS 



(10) 

(11) 
(12) 

(13) 



dV 



which allows us to obtain a conservative scheme by com- 
puting the fluxes of the advected quantity (p i; p n , PiVi, 



p n v n respectively in the equations above) at every inter- 
face on the grid and using the same flux to update ad- 
jacent zones. To solve the problem of calculating values 
of the advected quantities on the grid interfaces main- 



taining numerical stability (Stone & Norman 1992a), we 
used the second order van Leer interpolation method 



(van Leer 1977). This method is fast and accurate enough 
for the requirements of our simulations. As in ZEUS-2D 
the two-dimensional adve ction problem is simplified by us- 
ing directional splitting (Strang 1968), that means using 
two one-dimensional advection steps to construct the full 
solution. 

The boundary conditions used to solve the MHD equa- 
tions are determined by the physics of the phenomenon 
we are studying. We use periodic boundary conditions in 
y (note that with the definition of ip in equation (|^), the 
periodicity of ip ensures conservation of the total vertical 
magnetic flux through the simulation zone). At z — we 
launch an Alfven wave by giving the whole fluid (neutrals, 
ions and magnetic field lines) a motion in y which, in this 
first numerical test of the mechanism, will be limited to a 
single periodic oscillation: 

v y (y, z = 0) = v t cos(wi) . 

This will thus propagate upward as an Alfven wave. At the 
opposite boundary (z max ) we impose reflective boundary 
conditions; they have no real effect since the total length 
in z is chosen such that the waves are heavily damped 
before they reach that point. At the moment, in order to 
isolate the effect we want to study, we impose that there 
is no flux of matter from z = 0. In this manner we ensure 
that the variations in the ion or neutral densities do not 
result from inflow of matter at the lower boundary. On 
the other hand this results in severe constraints on the 
code because the wave pressure pushes matter upward, 
resulting in very low densities at the first grid point. This 
turns out to be the most stringent limit on the parameters 
we can use {e.g. the perturbed velocity Vt). 

The simulations have other numerical constraints. To 
ensure that our calculations are fully consistent the dimen- 
sions of our spatial grid must be large enough to allow the 
ion-neutral interaction, that is, larger than the ion-neutral 
(neutral- ion) mean free path, I = v/is, where v is the col- 
lision frequency. For the characteristic values used in the 
simulations the grid is more than 10 times the maximum 
mean free path for the neutrals, which is enough for the 
necessary ion-neutral interaction. 

We also wish to launch the wave in a manner such that 
it is damped not right in its emission zone (near z = 0, 
where we create it) but away from it, so that we can clearly 
identify the processes obtained. We do this by imposing a 
strong vertical density gradient. Thus the density at z = 
can be taken such that the wave initially propagates well 
(lu <C v n i — PPi), but later reaches an altitude where 
^ ~ v n i so that it is strongly damped and ambipolar pro- 
cesses can act. At higher altitude we maintain a small 
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Fig. 2. Left: Initial 2D distribution of the ionization fraction Z = pij p n in the calculation grid. Right: Ionization 
fraction after 6.7 x 10 3 Alfven times. The contrast along y is strongly enhanced at the altitude where the wave is 
damped. 



but constant density. Finally, in order to initiate the fil- 
amentation process we need a transverse density profile: 
the neutral density is initially independent of y, but the 
ion density has a horizontal profile, so that the ionization 
fraction is higher along the central field line (y — 0). The 
magnetic field is chosen to ensure MHD equilibrium in the 
horizontal direction, that is: 



d_ 

dy 



Pi 



8tt 







(14) 



Thus we expect that ions in the central flux tube will be 
compressed, and the neutrals expelled, by the ambipolar 
filamentation process. 

On the other hand, in order to save computation time 
we must take the largest time step that guarantees the 
numerical stability for the set of difference equations. Our 
time explicit code has to satisfy a Courant-Friedrichs- 
Lewy (CFL) stability condition, 



At < 



(15) 



derived applying a von Neumann stability analysis 
(Richtmayer & Morton 1957). Physically, this condition 
sets the higher limit of the distance that information 
can travel in one time step (waves or fluid motion) 
to be the minimum size of the discrete elements or 
'zones' in the spatial grid. For multidimensional systems, 
a suitable time step limit is the smallest of all one di- 
mensional CFL conditions in each coordinate direction 
( fBtone fc Norman 1992a| ). This gives strong constraints 



since (a) we need to consider the highest Alfven velocity, 
obtained where the density is lowest (b) the CFL condi- 
tion involves waves with the shortest wavelengths, i.e. high 



frequencies. At high frequency, as explained in the intro- 
duction, the relevant Alfven velocity involves only the ion 
density. Therefore compared with the Alfven velocity at 
z = the one used for the CFL criterion is higher by the 
density ratio (pi,max/ Pi.-min) 1 ' 2 along the vertical density 
profile, and by the inverse of the ionization ratio p n /pi- 
These two quantities are large and force us to use very 
short time steps, typically less than 10~ 2 of the Alfven 
time through a grid cell at z = 0. 

Let us now consider the time scales associated with the 
mechanism of ambipolar filamentation. Three time-scales 
appear. The first one is the period of the wave. The second 
one involves the response of the fluid to the wave pressure; 
in the basic mechanism of ambipolar filamentation, the 
neutrals feel a turbulent pressure: 



PTurb = ^Pn {vl) 



(16) 



(where the angular brackets mean a time average), which 
has a gradient along y associated with the gradient of pi\ 
they respond by a pressure perturbation to re-establish 
equilibrium. This is established after a time ~ L/c s , where 
L is the scale of the horizontal ionization inhomogeneity. 

However, at this stage the situation keeps evolving 
since the new neutral and ions density profiles are more 
peaked, resulting in a further peaking of PTurb- Thus, on 
a much longer time scale, the ionization contrast will keep 
growing dTagger et al. 1995[ ); in our simulations we find 
that typically ~ 10 2 Alfven times are needed to reach an 
equilibrium. Such long time scales mean that any small 
numerical resistivity will result in a diffusion of the mag- 
netic flux, so that the filamentation of the magnetic field 
associated with that of the ions cannot be observed. 
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Fig. 3. Top panels show the ionization fraction (left) ions density (center) and neutrals density (right) profiles in 
y direction at the altitude where the ionization fraction Z reaches its maximum. Bottom panels are the respective 
vertical profiles at the same point. Dotted lines show the initial values. 



3. Results of the Numerical Experiments 

The simulations were carried out using the numerical 
model described in the previous section. They follow the 
evolution of a 2-D (y-horizontal and z-vertical directions), 
two fluid system (ions and neutrals) with low ionization 
fraction. The gas is initially threaded by a vertical con- 
stant magnetic field and perturbed by horizontal waves ex- 
cited at the footpoints of the magnetic field lines (z = 0). 
Since no approximations were made except the invariance 
in x, the following parameters have to be imposed on each 
simulation: Bq, the initial magnetic field, c s , the sound ve- 
locity, p, the ion-neutral collisional coefficient, v t and to, 
the amplitude and frequency of the wave and the profiles 
of pi and p n in the y and z directions. 

In order to clarify the results presented in all figures, 
we must briefly explain the scaling and units used by our 
numerical code. We have normalized the parameters to 
the characteristic scales of the problem we are solving. 
Therefore we have taken as units the initial Alfven velocity 
at z — 0, va, and the Alfven time defined as the time 
that takes such an Alfven wave in crossing one wavelength 
Xa- 

In the simulation presented in Figs. 1-6 we used a 
61x200 numerical grid. It is longer in the vertical direc- 
tion in order to allow the complete damping of the wave. 
As it was established in the previous section, we assume 
initial equilibrium in a fluid vertically stratified and sup- 
ported by gravity, with ions and neutrals densities decreas- 
ing sharply (by a factor ~ 10) with z (shown in Fig. 1). 



Initially the neutral density is independent of y, but the 
ion density shows a small enhancement in the y direction, 
along which lie the perturbed velocities associated with 
the wave. The magnetic field is adjusted so as to ensure 
MHD equilibrium in the y direction, given by (^). The 
resulting initial ionization fraction Zq is constant over z 
but shows a maximum at y = 0, as shown in Fig. 2 (left). 
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Fig. 4. Vertical profile of the transverse perturbed veloc- 
ities of ions (solid) and neutrals (dotted). The decoupling 
of the two fluids begins at the altitude where the wave is 
damped. 

A high amplitude Alfven wave is launched from z = 0, 
by making the whole fluid (neutrals, ions and magnetic 
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field lines) oscillate in y at a single frequency ui; this 
perturbation propagates upwards as an Alfven wave. Its 
frequency is chosen so that the wave propagates with- 
out damping at the lower part of the simulation grid 
(lj <C Vni), but is strongly damped at intermediate al- 
titudes, where pi (and thus v„i) was taken to decrease 
sharply. Therefore the filamentation will occur only at 
the intermediate altitude where the wave is damped 



(Tagger et al. 1995) but still retains a strong perturbed 



velocity. In the simulation presented in Figs. 2-4 and 6, 
the perturbed velocity of the wave at z — is of the order 
of the sound velocity (v t = 0.7c s ), which is taken to be a 
half of the Alfven velocity. 
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Fig. 5. Vertical profile of the ratio between the final and 
initial ions density at y — yz mam ■ 

Fig. 2 (right) shows the spatial profile of the ionization 
fraction at the end of the calculation. The initial contrast 
in the horizontal direction has grown and shrunk at a cer- 
tain height. At higher altitude small traces of the wave 
still appear, together with a reversed profile along y. An 
explanation to this behavior will be given later in this sec- 
tion. In the vertical profile of the ionization fraction (Fig. 
3b) we note that the strong peak in Fig. 2 (right) is lo- 
cated at z ~ 32, actually where the wave begins to be 
damped (see Fig. 4) but still keeps enough amplitude to 
make the neutral expulsion mechanism effective. 

Therefore we find a strong amplification (by a factor 
~ 2.5) of the contrast in Z, at the altitude where the 
wave is damped (Fig. 3a). At the same altitude the ion 
density contrast along y, 5 pi, grows and the more ion- 
ized region shrinks (Fig. 3c), as expected. Moreover the 
ion density decreases on the sides of the peak, suggesting 
ion motion towards the most ionized regions. In Fig. 5 we 
can clearly quantify the final increase in ion density. We 
also see in the simulation how the neutrals are expelled 
from the most ionized regions (Fig. 3e) , generating in the 
density profile a 'central' minimum of the same order as 
the increase achieved by the ions. In all horizontal profiles 
of Fig. 3, the enhanced minimum/maximum are not cen- 



tered at y = because of the lateral motion due to the 
wave. On Fig. 4 we compare the ion and neutral velocity. 
A phase lag (corresponding to the damping of the wave 
by ion/neutral friction) appears around z = 32. 

There is an additional effect acting to expel the neu- 
trals from the most ionized regions. Since along those field 
lines the wave is less damped, its perturbed velocity at a 
given height is larger. However this contributes with a 
v.Vv term on the ions as well as on the neutrals, so that 
it should not contribute to the growth of the ionization 
contrast. 

Besides the strong peak in the spatial distribution of 
the ionization fraction, there is another feature which is 
a direct output of the enhancement of the ion density in 
the central magnetic field lines. The condensation of the 
ions must be accompanied by a compression of the mag- 
netic field lines at the same point. Then magnetic tension 
causes the compression of the field lines upwards. But as 
noted before, numerical resistivity has allowed the field 
to diffuse almost totally. A small residual intensification 
at the central field lines is left but it is invisible at low 
altitudes, where the wave dominates the field dynamics. 
However, at the highest grid zones, where the wave is al- 
most completely damped, the increase in B can be barely 
detected. That residual enhancement in B causes the rise 
of the magnetic pressure, which results in an expulsion 
of ions from the central field lines at high altitudes. This 
causes a reversed ion density profile, as observed in Fig. 
2. 

In our simulations the filamentation process seems only 
limited by two facts. Firstly, pi can reach zero (so that the 
code crashes) at large y. In this case the process is so effi- 
cient that all the ions at large y are evacuated towards the 
most ionized regions, at the altitude where the filamenta- 
tion occurs. Secondly, pi can reach zero in the lowest grid 
zones because the matter is pushed upwards by the wave 
pressure. Both processes impose severe restrictions to the 
model parameters, in particular to the highest value of v t 
we can use. The second limitation could be easily over- 
come by changing the boundary conditions to the more 
realistic case of allowing flux of matter at z = 0. However, 
as discussed above, we prefer not to do it in these first 
simulations since it would make it difficult to distinguish 
the enhancements of ion density due to the filamentation 
process from the ones due to this vertical flow. This will 
be necessary, on the other hand, in future realistic simu- 
lations of spicules. 

Another unwanted effect is visible in our plots, for the 
largest values of the turbulent velocities: perturbed quan- 
tities show a small-scale sawtooth appearance, of numeri- 
cal origin. They are associated with sharp features in the 
vertical velocity, although this component of the velocity 
remains much smaller than the horizontal one. This phe- 
nomenon is generated at the first vertical grid points, and 
results from a mismatch between the initial condition we 
impose (e.g. equal ion and neutral velocities at z = 0) 
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and the wave properties. This effect remains small how- 
ever and we have checked, by varying the grid size, that 
it does not alter our results. 

Fig. 6 shows the evolution of the maximum value of 
the ionization fraction for several initial wave amplitudes. 
For the lowest values of the perturbed velocity a state 
of equilibrium is achieved early in the calculation. How- 
ever, for the highest velocities the hlamentation is more 
efficient and causes the simulation to crash. In the sim- 
ulation shown in Fig. 3 (v t = 0.35va), this is due to the 
complete depletion of ions at the first grid zones, although 
very small values of pi are also achieved at high y, at the 
altitude where Z is maximum. In fact the curves corre- 
sponding to the highest velocities show an oscillating be- 
havior over t ~ lO 3 ^, which can be a consequence of the 
extremely low values of the ion density in those regions. 

Fig. 6 shows that the process acts faster for higher ve- 
locities. The asymptotic value of Z max varies quadratically 
with vt , as expected from the theory ( [Tagger et al. 1995| ) . 
However we find that this value (and accordingly the hor- 
izontal pressure gradients obtained for pi and p n ) depends 
on the initial profiles, i.e., more peaked initial profiles re- 
sult in more peaked final ones. Theory would lead us to 
expect that the final equilibrium, balancing the pondero- 
motive force of the wave with the pressure gradient of the 
neutrals, should be independent of the initial state. We 
believe that this may be due to the vertical transport of 
ions and neutrals from the lower and higher regions of the 
simulation grid. In these regions the filamentation process 
does not act, so that they retain a memory of the initial 
conditions and can feed the region of wave absorption with 
additional ions and neutrals. However we have been un- 
able to prove it with the present limitations, due to the 
boundary condition at z = 0. We thus defer the treatment 
of this question to future work where the detailed physics 
will be considered in more realistic conditions. 

4. Discussion 

The aim of the simplified 2-D simulations presented in this 
paper was to provide a first numerical test of the ambipo- 
lar filamentation mechanism. We have tried to limit the 
physics involved to the minimal ingredients needed, in or- 
der to clearly separate the expected physical effect (the fil- 
amentation) from other sources of variation of ion and/or 
neutral density. In particular we have limited ourselves to 
the excitation of a monochromatic wave and forbidden any 
inflow of matter from the lowest grid boundary, although 
this turned out to cause the most severe numerical limita- 
tion, by creating a vanishing ion density on the first grid 
points. 

We have found that the mechanism is strong and effi- 
cient, even with these constraints and simplifying assump- 
tions. In the most efficient of the simulations, the per- 
turbed velocity was close to the sound velocity taken to be 
a half the Alfven velocity (so the conditions are similar to 
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Fig. 6. Time evolution of the value of maximum ioniza- 
tion fraction for different initial perturbed velocities. 



that of the interstellar medium). However such velocities 
are still lower than the observed supersonic motions which 
would make the mechanism even more efficient, since we 
have found that its efficiency goes with the square of the 
perturbed ion velocities. We have obtained similar results 
in cases closer to the conditions of the solar spicules, where 
the sound velocity is a lower fraction of the Alfven velocity 
than in the results presented here (cs/va = -5). 

Future work will go in two directions: we will excite a 
more complete spectrum of waves from the lowest bound- 
ary. This should cause the filamentation to be less local- 
ized, since waves of different frequencies will be damped 
(and cause filamentation) at different heights. 

A more important evolution will be to relax the condi- 
tion of no vertical flow from z = 0. This will make the sim- 
ulation more efficient since we found that this boundary 
condition, chosen to clarify the physics, caused the most 
severe numerical limitation. It will also allow us to make 
more realistic simulations of solar spicules. In that case, 
we expect our mechanism to act together with the verti- 



cal flow of matter discussed in Haerendel (1995), Tagger 
et al (1995) and Dc Pontieu & Haerendel ( pL998| ). Upward 
traveling Alfven waves generated in the photosphere are 
expected to cause both ambipolar filamentation of the field 
lines and vertical acceleration of the gas, so that a bound- 
ary condition allowing matter to flow vertically from z = 
is necessary. First tests actually show an intense growth 
of the ionization fraction in the most ionized flux tubes. 
We also expect that, with these boundary conditions, we 
will be able to use higher perturbed velocities making the 
mechanism much more efficient. 

The effects of a more realistic equation of state and, 
of course, the inclusion of ionization and recombination, 
are also obvious extensions to this work. This should 
make us able to address more realistically the role of 
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our mechanism in the filamentary structure of the in- 
terstellar medium and the star formation process. In 
particular, the simulations shown in this work corre- 
spond to ionization fractions similar to those observed in 
the warm neutral component of the Interstellar medium 



( [Kulkarni fc Heiles 1987[ ). 

Models of star fo rmation (since the early work 
of Arons fe Max 1975| ) invoke molecular clouds sup- 
ported essentially by turbulent magnetic pressure and 
an ambipolar flow of the neutrals toward dense cores 
flArons fc Max 1975|; [gnu et al. 1987| |McKee et al. 1993 ; 
Vazquez- Semadeni et al. 199E; and references therein). In 
the molecular clouds conditions, where the turbulence is 
supersonic, we expect our mechanism to result in an ad- 
ditional pressure from the most ionized to the less ionized 
regions that efficiently separates the ions from the neu- 
trals, favoring the gravitational collapse of the latter. 

Another subject for future work could be the strong 
spatial and temporal intermittency of MHD turbulence. 
This is known t o occur in in tense current and vor- 
ticity sheets (s ee Spanglcr 1999|, and referenc es therein; 
Falgarone 1999 and Falgaronc fc Phillips 1990, for obser- 



vations in the interstellar medium and molecular clouds). 
We can expect our mechanism to be particularly efficient 
in these sheets where gradients of all quantities become 
strong on extremely small scales, and where turbulent mo- 
tions concentrate. 

The turbulent magnetic fields can be generated as a 
consequence of the Jeans-Parker instability during the 
molecular cloud formation process, or (as suggested by 
recent models) by the outflow lobes associated with 
young stellar objects; in the former case, large scale 
waves could feed the turbulent cascade in the cloud 
( |G6mez de Castro fc Pudritz 1992 ) . In the latter case, the 
outflows act as the sources of turbulence and the en- 
ergy released heats the cloud by the ambipolar drift 
( Nomura et al. 1999 ). In those contexts, there would be 
localized sources for the turbulence, that would propa- 
gate in a stratified medium, so that the physics studied 
here could be readily applied. On the other hand for a de- 
tailed study of the turbulent cascade and its effects on the 
interstellar medium, a pseudo-spectral technique would be 
better adapted to allow a random, non-localized excitation 
of the waves. The relevance of the role of the ambipolar fil- 
amentation in those conditions will be the object of future 
studies. 
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